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Abstract 

We demonstrate the existence of stable three-dimensional spatiotemporal solitons (STSs) in 
media with a nonlocal cubic nonlinearity. Fundamental (nonspinning) STSs forming one-parameter 
families are stable if their propagation constant exceeds a certain critical value, that is inversely 
proportional to the range of nonlocality of nonlinear response. All spinning three-dimensional STSs 
are found to be unstable. 

PACS numbers: 42.65.Tg, 42.65.Sf,42.70Df 



1 



Spatiotemporal solitons (STSs), also referred to as "light bullets" have attracted a 
great deal of attention in optics, see a recent review Q|. These are multidimensional pulses, 
which maintain their shape in the longitudinal (temporal) and transverse (spatial) directions 
due to the balance between the group-velocity dispersion, diffraction, and nonlinear self- 
phase modulation. However, solitons in media with the cubic self-focusing nonlinearity, 
obeying the nonlinear Schro dinger (NLS) equation, are unstable in two and three dimensions 
(2D and 3D), because of the occurrence of collapse in the same model (jj . Several possibilities 
to arrest the collapse were considered, such as periodic alternation of focusing and defocusing 
layers and various generalizations of this setting j^j], and the use of media exhibiting 
saturable [6 J or quadratic (x^) 7} nonlinearities; tandem structures, composed of alternating 
linear and x^ layers, were experiment in this field was 

the creation of quasi-(2+l)-dimensional STSs in bulk x^ samples 0,0]. Other theoretically 
developed approaches use off-resonance 
media jll]. 

Collapse does not occur either in x^ media whose nonlinearity is nonlocal Q|, therefore 
they may also give rise to stable solitons, see review [ljj. 2D spatial solitons stabilized 
by the nonlocality were observed in vapors ^| and lead glasses featuring strong thermal 
nonlinearity |l5||; in the latter case, elliptic and vortex-ring solitons were reported. Optical 
ID solitons supported by a nonlocal x nonlinearity were also created in liquid crystals 
Q. Further, self-focusing in photorefractive m edia Q, periodic lattices Q, vortices Q, 
and spatial solitons in soft matter [20] were considered in the context of nonlocality. 

The objective of this work is to demonstrate that one-parameter families of stable (3+1)- 
dimensional STSs are possible in media with nonlocal x nonlinearity. For this purpose, 
we consider a 3D model based on a general system of coupled equations for the complex 



field amplitude q and nonlinear correction to the refractive index n (see, e.g., 2l(); in a 
normalized form, the equations are 

iqt. + (l/2)(q m + g c? + Dq TT ) + qn = 0, 

d(n vv + n^) - n + \q\ 2 = 0. (1) 

Here, r], ( and £ are the transverse and longitudinal coordinates, scaled, respectively, to 
the beam's width and its diffraction length, r is the reduced temporal variable, and D is 
the ratio of the diffraction and dispersion lengths. We consider the case of anomalous 



temporal dispersion, D > 0, and then set D = 1 by means of an obvious scaling. Lastly, 
yd determines the correlation length A corr of the nonlocal nonlinear response [note that by 
setting d = one turns Eqs. (0) into the ordinary 3D NLS equation with self-focusing]. In 
fact, after setting D — 1, additional rescaling of Eqs. (0) makes it possible to set d — 1, so 
as to cast the system into the parameter-free form. Nevertheless, we prefer to keep d as an 
explicit parameter, as it directly controls the system's nonlocality degree. 

The nonlocal nonlinearity in Eqs. is typical for light propagation in liquid crystals and 
for thermal nonlinearity in optical media jla. For the derivation of the model equations 
((TJ), the usual approximation of the slowly varying amplitude is adopted, along with an 
assumption of fast temporal relaxation of the refractive-index perturbations, therefore the 
second equation in does not contain the term n TT . The latter assumption, which is 
essential for the physical justification of the nonlocal model including the temporal-dispersion 
term, implies that the relaxation time, r re i ~ A corr /c (c is the light speed), must be of the 
order of the pulse duration T of the STS to be constructed. Whether such conditions can 
be met in available optical media featuring nonlocal nonlinearities, such as liquid crystals 
which tend to exhibit slower relaxation times, is a question that remains a challenge. 

Equations (P) conserve the energy E — J J J \q(i],(, T )\ 2 drjd^dr, Hamiltonian H, the 
momentum P v ^ in the transverse plane, and angular momentum along the longitudinal 
direction. Stationary solutions to Eqs. (JTJ) are looked for as q = w(r, r) exp[i(&£ + S9)], 
n = n(r, r), where r and 9 are the polar coordinates in the (r], () plane, b is a real propagation 
constant, integer S is the vorticity ("spin"), and real functions w and n obey the equations 

(w rr + r^Wr + w TT ) - (2b + r~ 2 S 2 ) w + 2wu = 0, (2) 
d (n rr + r _1 n r ) — n + w 2 = (3) 

(for this solution, the angular momentum is proportional to the energy, M% = SE). 

We have numerically found families of localized solutions to these equations, dealing 
with the corresponding two-point boundary-value problem by dint of the standard band- 
matrix algorithm. Typically, grids with 241 x 240 and 201 x 360 points were used for 
the computations of the fundamental (S = 0) and spinning solitons, respectively. In Fig. 
1, we display the energy characteristic, E = E(b), for one-parameter families of the thus 
constructed soliton solutions. In view of the above-mentioned possibility to eliminate d by 
means of rescaling, panel (a) in Fig. 1 is, as a matter of fact, a blowup of a part of panel 
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(b) corresponding to < b < 2.4. A noteworthy feature is that the 3D solitons exist only 
above a finite energy threshold. 

The stabilizing effect of the nonlocality for the fundamental solitons is seen in Fig. 1: 
except for a narrow interval of small wavenumbers, < b < b CT , the solitons are expected 
to be stable, as they satisfy the Vakhitov-Kolokolov (VK) criterion, dE/db > 0, which is 
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a necessary (but, generally, not sufficient) stability condition for the soliton family 
(note that the instability of the 3D solitons in the local NLS equation precisely follows this 
criterion). The bending of the curve b = b(E) for fundamental 3D solitons, from the negative 
to positive slope, was also found in a different model fa medium with a Gaussian nonlocal 
kernel) by means of a variational approximation 23]. Below, it will be demonstrated that 
the VK criterion is actually sufficient for the stability of the fundamental solitons, but not 
for spinning ones, which are always unstable. 

For d — 1, the critical wavenumber which separates the negative- and positive-slope 
regions for the fundamental solitons in Fig. 1 is b^ ^ = 0.565, and, in view of the scaling 
invariance of the model, bir~°\d) = 0.565/ci. The energy of the 5 = soliton at the critical 
point is Ect~°^ = 42.60 for d = 1 (in the model with d ^ 1, E cr scales as Vd, see Fig. 
1). It also follows from here that Ecr = 32.02 ( bc T J , which may be compared to 
the scaling law for unstable fundamental solitons in the 3D local NLS equation, E^ S=QS) = 
C b-^ 2 , with C ^6.67 0. 

Figures 2(a,b) display the shapes of typical stable and unstable fundamental soliton for 
d = 10, at a fixed value of the energy (E = 140). It is seen that high-amplitude solitons are 
stable, whereas low-amplitude ones are unstable. We notice that both the low- and high- 
amplitude solitons feature spatiotemporal ellipticity, being broader in space than in time, 
which can be easily explained by a perturbative analysis of Eqs. for small d. Typical 
shapes of unstable S = 1 solitons are shown in Figs. 2(c,d). 

Full stability of solitons was investigated using the equations for small perturbations 
linearized around the stationary solution. Accordingly, solutions including perturbations 
with an infinitesimal amplitude e are looked for as 

q = e ib ^ +lSe {w(r, r) + e [f(r, T)e^ +ije + g*(r, r)e 5 ^- ije ] } , (4) 
n = n {r, r) + e [p(r, r)e^ +ije + p*(r, T)e s ^~ lje ] , (5) 

where J is an arbitrary integer azimuthal index of the perturbation, 5 is the instability 



growth rate, the asterisk stands for the complex conjugation, and the eigenfunctions /, g 
and p obey the equations 

2i5f + f rr + r~ l f r + f TT - [2b +(S + Jfr- 2 ] f + 2(nf + wp) = 0, 
-2i5g + g rr + r~ l g r + g TT - [2b + (S - J) 2 r~ 2 ] g + 2{ng + wp) = 0, (6) 
d(p rr + r~ x p r ) - (1 + d,J 2 r~ 2 )p + w(f + g) = 0. 

The growth rate 5 was found as an eigenvalue at which Eqs. (JHJ) has a nonsingular 
localized solution. In Figs. 3(a), 3(b,c), and 3(d), we plot the instability gain, Re(<5), vs. 
the propagation constant b of the unperturbed solution, for the STSs with S = , 1, and 
2, respectively. The stable solitons are those for which Re(5) = for all (integer) values 
of J. Due to the scaling invariance of Eqs. (JTjl, the curves in a given panel, pertaining to 
different values of d, are obtained from each other by the scaling transformation, therefore 
they actually display essentially the same stability and instability regions, but on different 
scales. Figure 3(a) shows Re(<5) for J = 0, as J ^ does not yield any instability for the 
fundamental solitons; the stability region revealed by this figure, b > 6 cr , is precisely the 
same as predicted above by the VK criterion. For the spinning solitons, the perturbations 
with J > 1 destabilize a part of the families, see, e.g., Fig. 3(c) for S = 1 and J = 2 , but 
entire STS families with S > 1 are unstable against the perturbations with J = 1, see Fig. 
3(b,d). The latter instability mode implies a trend to splitting of the vortex soliton into a 
set of two fundamental ones |2|, which is corroborated by direct simulations below. 

Note that stable 3D spinning optical solitons (with S = 1) were only found in media with 
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competing focusing and defocusing nonlinearities, viz., '■ or '■ X^ 
other hand, stable 2D spatial (rather than spatiotemporal) vortex r ings were experimentally 
observed in a medium featuring the thermal nonlocal nonlinearity 

The predictions of the linear stability analysis were checked in direct simulations of Eqs. 
(0), which were run by means of a standard Crank- Nicholson scheme. The nonlinear finite- 
difference equations were solved using the Picard iteration method, and the resulting linear 
system was handled with the help of the Gauss-Seidel iterative procedure. To achieve good 
convergence, we needed typically, twelve Picard's and four Gauss-Seidel iterations. The 
initial conditions for perturbed solitons were taken as g(£ = 0) = w(r), (,t)(1 + ep), and 
n (£ = 0) = 71(77, C, t)(1 + ep), where e is, as above, a small perturbation amplitude, and 
p was either a random variable uniformly distributed in the interval [—0.5,0.5], or simply 



p = 1 (uniform perturbation). 

First, we have checked that all the fundamental STSs that were predicted above to be 
stable (stable branches are shown with solid lines in Fig. 1), are stable indeed against random 
perturbations; Fig. 4 displays an example of self-healing of a stable soliton with the initial 
perturbation amplitude e = 0.1. A small uniform perturbation (p = 1) applied to a stable 
soliton excites its persistent oscillations, which suggests the existence of a stable intrinsic 
mode in the soliton. On the other hand, direct simulations show that those fundamental 
solitons that were predicted to be unstable decay into radiation, if slightly perturbed. We 
also simulated self-trapping of a stable fundamental soliton from an initial spatiotemporal 
pulse of an arbitrary form. An example is shown in Fig. 5 for an isotropic Gaussian input, 
which generates an anisotropic (elliptic) soliton. We also simulated the evolution of unstable 
spinning solitons. Most typically, they follow the prediction of the linear stability analysis 
and split into two stable fundamental solitons, see an example (for S = 1) in Fig. 6. 

In conclusion, we have demonstrated that nonlocal cubic nonlinearity is sufficient to sta- 
bilize 3D solitons, which suggests a new approach to making of 3D spatiotemporal solitons 
in optics, which thus far evaded experimental observation. The stability of the funda- 
mental solitons was demonstrated through the computation of the corresponding stability 
eigenvalues, and in direct simulations. Their robustness and, hence, physical relevance was 
demonstrated by self-trapping from arbitrary input pulses. On the other hand, all the spin- 
ning 3D solitons in the nonlocal medium are unstable against splitting into a set of stable 
fundamental solitons. 

This work was partially supported by the Government of Spain through grant BFM 2002- 
2861, by the Ramon- y-Cajal program, and by Deutsche Forschungsgemeinschaft (DFG), 
Bonn. The work of B.A.M. was supported, in a part, by the Israel Science Foundation 
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FIG. 1: Energy E vs. propagation constant b for soliton families with different vorticities, 5 = 0, 1, 
and 2 (numbers labeling the curves), at different values of the range of nonlocality of nonlinear 
response \fd: (a) d = 1, (b) d = 10. Full and dashed lines show stable and unstable solitons. 

FIG. 2: Cross-section shapes of typical 5 = and 5=1 solitons in the transverse (r) and temporal 
(r) directions, (a) and (b): 5 = 0, d = 10, E = 140; full and dashed lines correspond to stable 
and unstable solitons, with b = 0.15 and b = 0.021, respectively, (c) and (d): 5 = 1, b = 1; solid 
and dashed lines correspond to d = 1 and d = 100, respectively (due to the scaling invariance, the 
latter is tantamount to d = 1 and b = 100). 

FIG. 4: Isosurface plots illustrating the stability of a fundamental soliton corresponding to d = 10, 
6=1, and E = 178. (a) and (c): the initially perturbed soliton, at £ = 0; (b) and (d): the 
self-cleaned one at £ = 360. Here and in Figs. 5 and 6, the upper and lower rows show \q\ 2 and n, 
respectively. 



FIG. 5: Self-trapping of a fundamental soliton, for d = 1: (a) and (c) - the initial Gaussian pulse 
with energy Eq = 54; (b) and (d) - the soliton at £ = 60. 



FIG. 6: Splitting of an unstable 5 = 1 soliton with d = 100 and b = 0.05. (a) and (d) £ = 0, (b) 
and (e) £ = 1400, and (c) and (f) £ = 1600. 



FIG. 3: The real part of the perturbation growth rate vs. the propagation constant of the unper- 
turbed soliton. (a) 5 = 0, J = 0, for d = 1 and d = 10, (b) 5 = 1, J = 1, (c) 5 = 1 , J = 2, (d) 
5 = 2, d = 20. In (b) and (c), values of d, and in (d), values of J are indicated near the curves. 
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